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Abstract In this paper we study the effect of the anisotropic stress gen- 
erated by neutrinos on the propagation of primordial cosmological gravita- 
tional waves. The presence of anisotropic stress, like the one generated by free- 
streaming neutrinos, partially absorbs the gravitational waves (GWs) propa- 
gating across the Universe. We find that in the standard case of three neutrino 
families, 22% of the intensity of the wave is absorbed, in fair agreement with 
previous studies. We have also calculated the maximum possible amount of 
damping, corresponding to the case of a flat Universe completely dominated 
by ultrarelativistic collisionless particles. In this case 43% of the intensity of the 
wave is absorbed. Finally, we have taken into account the effect of collisions, 
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using a simple form for the collision term parameterized by the mean time 
between interactions, that allows to go smoothly from the case of a tigthly- 
coupled fluid to that of a collisionless gas. The dependence of the absorption on 
the neutrino energy density and on the effectiveness of the interactions opens 
the interesting possibility of observing spectral features related to particular 
events in the thermal history of the Universe, like neutrino decoupling and 
electron-positron annihilation, both occurring at T ~ 1 MeV. GWs entering 
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the horizon at that time will have today a frequency v ~ 10 9 Hz, a region 
that is going to be probed by Pulsar Timing Arrays. 
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1 Introduction 

The presence in the Universe today of a stochastic background of gravitational 
waves (GWs) is a quite general prediction of several early cosmology scenarios. 
In fact, the production of gravitational waves is the outcome of many processes 
that could have occurred in the early phases of the cosmological evolution. 
Notable examples of this kind of processes include the amplification of vacuum 
fluctuations in inflationary [1] and pre-big-bang cosmology scenarios [2] , phase 
transitions [3], and finally the oscillation of cosmic strings loops [4]. In most 
of these cases, the predicted spectrum of gravitational waves extends over a 
very large range of frequencies; for example, inflationary expansion produces 
a flat spectrum that spans more than 20 orders of magnitude in frequency, 
going from 10~ 18 to 10 9 Hz. 

The detection of such primordial gravitational waves, produced in the 
early Universe, would be a major breakthrough in cosmology and high energy 
physics. This is because gravitational waves decouple from the cosmological 
plasma at very early times, when the temperature of the Universe is of the 
order of the Planck energy. In this way, relic gravitational waves provide us 
a "snapshot" of the Universe near the Planck time, in a similar way as the 
cosmic microwave background (CMB) radiation images the Universe at the 
time of recombination. 

The extremely low frequency region {y§ < 10~ 15 Hz) in the spectrum of 
primordial gravitational waves can be probed through the anisotropics of the 
CMB. In particular, gravitational waves leave a distinct imprint in the so-called 
magnetic or B-modes of its polarization field [5,6]. The amplitude of the pri- 
mordial spectrum of gravitational waves is usually parameterized through the 
tensor-to-scalar ratio r, i.e., the ratio between the amplitudes of the initial 
spectra of the tensor and scalar perturbations in the metric. The Planck satel- 
lite [7] is expected to be sensitive [5] to r > 0.05, correspondingto a density 
parameter J? gw (^) = (1/ p c )dp gw /d\ogv as faint as ~ 3 x 10~ 16 /i~ 2 (h is the 
dimensionless Hubble constant); future experiments are expected to enhance 
the sensitivity of two orders of magnitude [S]. 

On the other hand the operating large-scale interferometric GW detectors, 
although designed with the aim to detect astrophysical signals, can possibly 
also detect signals of cosmological origin [9 . They give complementary infor- 
mation with respect to the CMB polarization field since they probe a different 
region in the frequency domain. In particular the ground-based interferome- 
ters, such as the LIGO [TU], VIRGO [II], GEO600 [H] and TAMA300 [B] 
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experiments, operate in the range between 10 Hz and few kHz, and are ex- 
pected to be sensitive to fi gvJ h 2 > 10~ 2 . Even more interesting is the LISA 
space interferometer [14], that will hopefully operate in the 2020s. Not being 
hampered by the Earth seismic noise, it will probe the frequency region be- 
tween 10~ 4 and 1 Hz and will in principle be able to detect f2 gw h 2 > 10 -12 
at vo — 10 -3 Hz. According to theoretical predictions, a large enough GW 
signal at this frequencies can be produced, with the appropriate choice of pa- 
rameters, by a pre-big-bang accelerated expansion, by the oscillation of cosmic 
strings, or by the electroweak phase transition occurring at T = 300 GeV. Fi- 
nally, pulsar observations can be used to obtain information on the stochastic 
GW background, through the technique known as pulsar timing. The so-called 
"Pulsar Timing Arrays" will probe the region around 10~ 9 Hz [15] . 

In order to compare the theoretical predictions with the expected instru- 
ment sensitivities, one needs to evolve the GWs from the time of their produc- 
tion to the present. It is often assumed that gravitational waves propagate in 
vacuum, i.e., they freely stream across the Universe. In this case, the only ef- 
fect on a propagating GW is a change in frequency (corresponding to the usual 
redshift of the wavelengths caused by the expansion of the Universe) , and a 
corresponding change in the energy of the wave. However, GWs are sourced by 
the anisotropic stress part of the energy- momentum tensor of matter, so that 
the vacuum approximation is well-motivated only when this can be neglected. 
It is already known that the anisotropic stress of free streaming neutrinos acts 
as an effective viscosity, absorbing gravitational waves in the low frequency 
region, thus resulting in a damping of the B-modes of CMB [Tfflll7lH8U19ll2"0"] . 

In the present work we aim to start a study concerning the possible effects of 
the presence of neutrinos in other frequency ranges, like those probed by Pulsar 
Timing Arrays or interferometers. In particular, we aim to understand if events 
occurring in the thermal history of the Universe, like neutrino decoupling or 
the electron-positron annihilation, can leave an imprint in the spectrum of 
cosmological gravitational waves. The rationale behind this is that these events 
give rise to sharp changes in the neutrino density and in the neutrino mean 
free path; corresponding in turn to sharp changes in the anisotropic stress of 
the cosmological fluid. This would point to the fact that GWs entering the 
horizon before or after these events would experience a different amount of 
absorption. 

The paper is organized as follows. In Sec. 2, we introduce the basic equa- 
tions. In particular we introduce the coupled Einstein-Boltzmann system and 
recast it in a form that is very suitable for numerical integration. In Sec. 3, 
we show the results of the numerical integration of the Einstein-Boltzmann 
system, showing the time dependence of the GW amplitude and computing 
the amount of absorption for different values of the neutrino density and of 
the neutrino mean free-path. Finally, in Sec. 4 we draw our conclusions and 
put forward some ideas for the future. 
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2 Basic equations 

We shall use, all throughout the paper, natural units in which c = H = fcg — 1- 
Let us consider a gravitational wave, propagating on the background of a 
flat Friedmann Universe. In synchronous gauge, the spatial components of the 
perturbed metric are written as 0: 

9ij = a 2 (i) + (1) 

while the other components are left unperturbed: goo = — 1 and goi — 0. We 
will consider only the transverse traceless part of , representing a GW. Here 
a(t) is the cosmological scale factor, that evolves according to the background 
Friedmann equation: 

da\ 2 8ttG , 

*) =— aP (2) 

where p is the background density of the cosmological fluid (in general, we will 
use overbars to denote background quantities). 

The components of the tensor hij evolve according to [21] : 

d t^ + ( ~) d t h tJ - = l^G^f , (3) 



lJ ^\adt) ut l ° \a* J 13 ~ a 2 

where 77^ is the anisotropic stress, i.e. the traceless part of the three dimen- 
sional energy-momentum tensor TJ- of the cosmological fluid. It is then defined 
through the relation 77] = Tj - <5jT*/3 = Tj - V5], where V is the total 
pressure of the fluid (including possibly a small perturbation with respect to 
the background). 

The macroscopic properties of the cosmological fluid can be derived by the 
phase space distribution of its particles. The phase space is described by three 
positions x % and by their three conjugate momenta Pi = mdxi/ds. The proper 
momentum pi — p l measured by an observer at a fixed spatial coordinate is 
related to Pi by Pi = a(Sij + \hij)p> . As usual, the phase-space distribution 
function (DF) f(x z , Pj,t) of the particles gives the number of particles inside 
the 6-dimensional volume element: 

f(x\ P J ,t)d 3 xd 3 P = dN. (4) 

The energy-momentum tensor can be written in terms of the distribution 
function as follows: 

If P^P 
TZ = —jf (x l ,P jt t) -^f dP x dP 2 dP 3 , (5) 

where g is the determinant of the metric. 

In the following, we will use the comoving three-momentum qi — api in 
place of Pi as a momentum variable, and write it as qi — qrii where rij is a 
unit vector. We will also use the conformal time ry, defined by dt = a dr\ as our 



1 We use the ( h H — (-) signature for the metric. 
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time variable. Then we write / = f(x l ,q,rij,r)). Finally, we write the DF as 
the sum of a zeroth-order, unperturbed part /o, and a small perturbation: 

f(x i ,q,n j ,r 1 ) = fo(q)[l + V(x i ,q,n j ,T)] . (6) 

Using the fact that y/=g = a 4 /(l - \h) and dP 1 dP 2 dP 3 = (1 + \h)q 2 dqdQ, 
where h is the trace of hij and df2 is the infinitesimal element of solid angle 
around n, we can write: 

Tj = a- 4 j ^nS/o(9)(l + V)q 2 dqdf2 (7) 



where e = ^/q 2 + a 2 m 2 . We warn the reader that, even if we follow the con- 
vention of distinguishing between covariant and contravariant indices, the fact 
that we use the flat metric 8ij to raise and lower the indices of pi (and hence 
rij) means that equalities like that in Eq. ((7]) are not covariant. The unper- 
turbed phase space distribution is given by a thermal equilibrium distribution, 
i.e. by a Fermi-Dirac or Bose-Einstein distribution: 

/ ° (g)= (2 g 7 rVe^±l ' (8) 

where g s is the number of quantum degrees of freedom, and To is the present 
temperature of the particles. 

The DF evolves according to the Boltzmann equation: 

Uf] = C[f] (9) 

where the L = Df /Drj is the Liouville operator, and C is the collision oper- 
ator accounting for collisions between particles. Using the geodesic equation, 
the Liouville operator L can be cast in the form (to first order in perturbed 
quantities): 

l [f] -Df _9f { dx* df 1 i nJ -dft»j df (1Q) 
Drj drj drj dx % 2 dr\ dq 

Once the collision term is also specified, Eqs. (|2|), ([3]), (J7]), (HJ) and (TlO|) are all 
that is needed, at least in principle, to follow the propagation of a GW. 



2.1 Multipole formalism 

In this subsection, we will rewrite the coupled Einstein-Boltzmann system 
derived above to a form that is more suitable for numerical integration. With 
very small variations, this is the same procedure used when dealing with scalar 
perturbations [2 23. First of all, we note that we will only be concerned with 
massless particles as a source for the anisotropic stress, so we will set the mass 
m equal to zero in all the formulas derived above. Although we will be referring 
to these particles as "neutrinos" , for the purpose of computing the evolution 
of cosmological GWs they could actually be everything as long as they are 
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effectively massless, i.e. as long as the temperature of the cosmological plasma 
is much larger than their rest mass. 

Firs of all, we Fourier transform the spatial dependence of all the relevant 
quantities introduced above. With a slight abuse of notation, we shall use 
the same symbol to denote a given quantity and its Fourier transform. The 
Boltzmann equation in fc-space reads (dots denote derivatives with respect to 
conformal time): 



V + ihn l V - iriV^^r^ = ( n ) 
2 d\nq f 



In the case of massless particles, the dependence of the DF from q can be 
integrated out. In particular, after defining 

FAk, n,, r) = HMf^l y.^ , (12) 

J q i fo{q)dq 

we can multiply the Boltzmann equation by q 3 fo and integrate over q\ the 
result is: 

F„ + i hn l F v + 2^nV - f q 3 C[f]dq, (13) 

where we have also used the fact that p v = 4na~ 4 J q 3 fo(q)dq. Then we define: 

Fij{ki, H,t) = J [niTij - -^-J F v dcj), (14) 

where cj> is the polar angle, so that the infinitesimal solid angle element dfi = 
sin 8d8d(f>. Multiplying Eq. (|T3|) by (mrij — <5y/3) and integrating over 0, we 
get (we define as usual fj, = k ■ h): 

Fij + ikfiFij + 2hi m J n l n m f n,^ ^-J dcj) = (15) 

where we have defined the "collision term" Cij : 

dj = ^=J q 3d( i J J *t> {^i n i - 5 -f) d ^ > ( 16 ) 

Now we expand Tij and Cij in Legendre polynomials: 

oo 

F l3 (h, /j,, t) = + T )P t {n) (17) 

£=0 
oo 

djih, Ml r) = ^Hj^^ + ljcg^fci, t)P^) (18) 

£=0 
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In order to obtain the "tower" of (infinite) differential equations for the T\ 
we multiply Eq. (fTS"|) by (i /2)Pi, integrate over /j, and use the orthogonality 
relation of the Legcndre polynomials, i.e: 



.! p 2i + l 



(19) 



The detailed calculation is shown in the appendix. The final result is: 



t-(o) 



T, 



,(0) 



(2) _ 



15 



3^? - 2j£> 



16tt • 
105 ij 



r (2) 
'°y ' 



(20) 
(21) 



5J*. 



(5) 



4J". 



(3) 



315 



,(4) 



(22) 



J 7 , 



« _ 



2^+1 



(£ + i)^ +1) -i?./; 



(<-i) 

y 



(^0,2,4). (23) 



This system of infinite first-order ordinary differential equations is completely 
equivalent to the original Boltzmann equation and can be solved with fairly 
standard numerical methods for ODEs. 

The system should be closed with the evolution equation for hij, namely 
Eq. In Fourier space, and using conformal time, this reads: 



h{j + 2T~Lhij 



k 2 h l3 = lGnGIIij, (24) 



where T-L is the "conformal" Hubble constant H — a /a. The anisotropic stress 
Hij is given by: 



- ^<°> (25) 



In deriving this expression we have used the fact that J niTijdfi = 4atSij/3. 
Then we finally get: 



hij + 2Hhij 



k 2 hij = 4Ga?p v r i3 



(0) 



(26) 



Finally, we have to specify the initial conditions for the integration. By 
studying the behaviour of the solution when the wave is far outside the horizon 
(krj > 1) it can be seen that the right initial conditions are hij — and 



w _ 



0. The initial value of hij is arbitrary but the equations can 



always be rescaled to have h 



(()) 
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3 Interaction of gravitational waves with neutrinos 

We can now use the equations derived in the previous section to study the 
effect of neutrinos on the propagation of cosmological GWs. We restrict our 
attention to waves entering the horizon well before the time of matter-radiation 
equality, corresponding to a redshift z ~ 10 4 . This corresponds to waves with 
a present frequency v 3> 10 -16 Hz. During the radiation-dominated era, a oc i] 
and H = l/rj. It can be shown that using the time variable u=krj the evolution 
equations can be recast in a form such that k does not explicitly appear, so 
that the evolution in with respect to u is independent from the frequency of 
the wave. Also, the convenience of using u is that u = kr/ <~ 1 corresponds to 
the time of horizon crossing. The results of the numerical integration should 
be compared with the solution in the absence of anisotropic stress (77^ = 0), 
i.e. hij = h\j smu/u. 

(i) 

First we consider the case of a vanishing collision term, C\j = 0. This 
is the case after neutrino decoupling, occurring when the temperature of the 
Universe is T ~ 1 MeV and 1 + z ~ 10 10 . At lower temperatures, neutrinos do 
not interact with the other particles in the cosmological plasma so that they 
are free-streaming and collisions are effectively absent. This is basically the 
case that was considered in Ref. |18) . In order to parameterize the neutrino 
density, it is useful to introduce the quantities R v = p u /p^ and f v = p v /p, 
related by f v = R v /(1 + R u ). Taking the standard case of three neutrino 
families with a temperature T v = (4/H) 1 ' 3 ^, we have that: 



so that R v = 0.6813 and f v — 0.4052. The results of the numerical integration 
performed using this value of f v is shown in Fig. [1] In the top panel we plot the 
evolution of hij (divided by its initial value h^) with respect to u (red curve). 
We see that, as it should be expected, the amplitude hij is constant outside 
the horizon and starts decreasing after the horizon crossing; this is mainly 
due to the redshift caused by the expansion of the Universe. However, when 
comparing with the zero-stress solution smu/u (black short-dashed curve) 
it can be noticed that the wave suffers an additional damping, caused by the 
anisotropic stress of the neutrinos. This is made even more clear in the bottom 
panel, where we plot the combination u 2 \hij\ 2 in order to see the behaviour of 
the intensity \hij\ 2 once the expansion of the Universe has been taken out. We 
se that the amount of damping with respect to the zero-stress case tends to a 
constant value. For /„ = 0.4, the intensity of the wave is 0.78 times its value 
in the absence of stress, so that roughly 22% of the intensity of the wave is 
absorbed. 



The amount of absorption increases with the neutrino fraction f v . It is worth 
stressing, at this point, that the cosmological neutrino background has not 




(27) 
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0.1 0.5 1.0 5.0 10.0 50.0 

u = kij 



Fig. 1 Top panel: Evolution of the normalized wave amplitude hij/h^j' for /„ = 0.4 (red 
solid curve), 1 (blue dashed curve), (black dotted curve). Bottom panel: The same as the 
top panel, but for the normalized wave intensity \hij/h^\ 2 times u 2 , so that the redshift 
due to the expansion has been taken out. 

been directly observed yet. Although strong deviations from the standard sce- 
nario with f v = 0.4 are unlikely, the possibility that /„ has a different value 
(due for example to the presence of additional particles in the early Universe) 
should be taken into account. For this reason we show in Fig. [T]also the ex- 
treme case f„ = l (blue long dashed curve), i.e. p v — p, meaning that all 
the matter content of the Universe is made of non-interacting, ultrarelativistic 
particles. This gives the maximum possible amount of damping: the intensity 
of the wave is nearly halved, being 0.57 times its zero-stress value. 

Then we turn to consider the effect of collisions. The exact computation of 
the collision term depends on the details of the interaction. However, a useful 
although rough approximation consists in writing C[f] = — fo&fr, where r is 
the average time between collisions. This will give c|f = —J^ jr. The key 
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parameter in defining the strength of the interactions is fcr, that is basically 
the ratio of the frequency of the wave to the frequency of the collisions. A 
very small value of fcr would correspond to very frequent and thus effective 
collisions; the right-hand sides of Eqs. (f2T) |) - (j!?5j) would be dominated by the 
collision terms and would have the solution oc e _I, / r , meaning that the 
anisotropic stress would decay exponentially. On the other hand, a large value 
of fcr would correspond to rare collisions and in the limit fcr — > oo the colli- 
sionless result should be recovered. Without resorting to any particular model, 
we show in Fig. [5] the result of the numerical integration for different constant 
values of fcr ranging from 0.5 to 10. We fix the neutrino fraction to the extreme 
value f v = 1 in order to make the differences more evident. We see from the 
figure that the smaller the value of fcr, the more the amplitude of the wave 
tends to its undamped value of 1, while when the collisions are rare (large fcr) 
we recover the large damping found above. 




u = kjj 

Fig. 2 Evolution of the wave for f n u = 1 and different value of the mean collision time r. 
From top to bottom: kr = 0.5, 1, 2, 10. 



4 Conclusion and prospects 

In this previous sections we have studied the effect of the anisotropic stress 
generated by massless particles ( "neutrinos" ) on the propagation of cosmolog- 
ical gravitational waves. In particular, we have put the relevant equations in 
a form that is very suitable for numerical integration and allows a quite clear 
quantitive understanding of the effect of the relevant parameters. The pres- 
ence of anisotropic stress, like the one generated by free-streaming neutrinos, 
partially absorbs the GWs propagating across the Universe. In the standard 
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case of three neutrino families, the wave is damped by a factor 0.78 in intensity 
(0.88 in amplitude). We have also calculated the maximum possible amount of 
damping, corresponding to the case of a flat Universe completely dominated 
by ultrarelativistic collisionless particles. In this case the wave is damped by a 
factor 0.57 in intensity (0.75 in amplitude). Finally, we have taken into account 
the effect of collisions, and, using a simple form for the collision term param- 
eterized by the mean time between collisions, we have shown that we can go 
smoothly from the case of a tigthly-coupled fluid to that of a collisionless fluid. 

The dependence of the amount of damping from the neutrino fraction and 
from the effectiveness of the interactions opens interesting possibility. In fact, 
neither one of these quantities is really constant during the history of the 
Universe. The density of ultrarelativistic particles experiences several abrupt 
changes during the cosmological evolution, corresponding to the creation of 
new particle species when the temperature of the Universe is large enough. 
The effectiveness of the interactions depends on the number density of tar- 
get particles and on the interaction cross section, and both these quantities 
are a function of the temperature. Since the evolution of a gravitational wave 
is mainly affected by events occurring around the time of its horizon entry, 
and this in turn is directly related to the wave frequency, this open the in- 
teresting possibility of observing spectral features related to particular events 
in the thermal history of the Universe. A very promising frequency range is 
the nanohertz range. This roughly corresponds to GWs entering the horizon 
when T ~ I MeV, corresponding not just to one, but to two notable events in 
the thermal history of the Universe: neutrino decoupling and electron-positron 
annihilation. The first corresponds to the transition, for the neutrinos, from 
being a tigthly coupled fluid to being a collisionless gas, so it can be roughly 
thought as representing a (quite fast) change in the value of r. The second 
event, e + e~ annihilation, occuring shortly thereafter, changes the ratio of the 
photon and neutrino temperatures, and thus marks a sudden change in the 
value of f v from 0.72 to 0.40. Interestingly, the nanoHertz frequency region 
is going to be probed by the so-called Pulsar Timing Arrays [15] and thus 
represents a very promising observational target possibly allowing to increase 
our knowledge of the thermal history of the early Universe. 



Appendix 

In this appendix we show the calculations leading to Eqs. (| 2 01) - ((23)) from Eq. 
(ITSl) . Let us consider separately the three terms in Eq. (fT5"|) . 

First term. This is simply: 




(28) 
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Second term. We need to calculate 



— ik J d/j./j.FijPe 



(29) 



Using the recurrence relation: 

(e + i)p t+1 -(2e + i)nPt + ep t - 1 = o (30) 

we can express iiPeXp) in terms of Pe+i and Pt-\. We then get: 



i e+1 k 



f 1 i i+1 k ^ , n Z" 1 

/ d[i[iTijPi = -— (-0 m (2m+l)J-g" ) / d^ 



2^+1 



p — 



m=0 



2(£+l) 



(2*+l)(2m+l) 



M-l,m + 



2£ 



(2*+l)(2m+l) 



5/- 



2^+1 



Hr^< 1) +Hr i < 1) 



2^+1 



(31) 



TTwrd £erm. We need to calculate: 



y <fyt 2h lm J 



2hi m I d(j>n l n m I r^n; — 



(32) 



Let us start from the angular integration. It can be shown that, if Aij is a 
generic (three-dimensional) symmetric, transverse traceless tensor, then: 



Ai m / d<pnin m = 



I 

L 



>nin m n t nj = -(^ 2 - l) 2 Aj 



(33) 
(34) 



So the term reduces to : 



TT / d M 



2/l/r, 



2tt 



p^/x) = — ^m^-i) 2 ^^) = — 

(35) 



where the 7^) are related to the coefficients of the expansion in Legendre 
polynomials of the function (/1 2 - l) 2 = (8/15)P - (16/21)P 2 + (8/35)P 4 and 
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are given by: 



M) - 



16 
15 

32 
105 

16 
315 



I = 0, 



(36) 







= 4, 



otherwise. 



Putting all together we finally get 

k 



21+1 



ij 



(37) 



Note added in the arXiv version After this paper was published, we became 
aware that the effect of sudden changes in the radiation energy density on 
the spectrum of GWs has been studied in Ref. [23]. Moreover, in Ref. [24] the 
effects of free-streaming neutrinos have been computed up to second order in 
perturbation theory. 
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